Final Project

Sustainable Soil Management Practices and Impacts

Abstract

Industrial agriculture is threatened by a changing climate. With the current methods under threat and the environmental factors becoming increasingly difficult to mitigate, solutions are needed. One area of agriculture that is taking heavy damage is soil health. With pre-existing farming methods like tilling, fertilizer, and pesticide use, soil health is in steady decline. In this study, we identify the croplands of Bismarck, North Dakota, and determine how sustainable agriculture methods affect soil quality. Using indicators like soil quality index scores (SQIs),  Beta-Glucosidase Activity (BGLU), and Bacteria Maker FAME (BACSUM), we utilized a machine learning workflow to determine how the practices affect soil health. Our results suggest that cover cropping, diversified rotations, and grazed continuous cropping support better soil health compared to fallow-based methods. Traditional rotation practices do not show an increase in soil recovery and restoration of organic matter. While these practices are not detrimental to existing agricultural land, the switch to more sustainable alternatives is imperative for adapting to the outcomes of climate change. Recognizing how soil health responds to different management practices can serve as a starting point for further research in different agrosystems outside of North Dakota.

Introduction/Hypothesis

As climate change worsens, attention to industrial agriculture increases. The specific methods of the industry, like tilling, monocropping, and fertilizer and pesticide use, pose a threat to soil health (Belete & Yadete (2023)). These practices, while they may be efficient and produce steady yields, have further environmental implications that are not sustainable in the long run (Worley (2024)). Soil health is critical to food and water security, as well as overall ecosystem function. In ecosystems, soil organic carbon (SOC) is fundamental in the processes of water permeation, nutrient retention, and aggregate stability. As Zhao et al. state, “Even a small change in the SOC can substantially affect not only the climate but also the stability of ecosystems, because of its decisive role in the exchange of carbon between the soil and atmosphere and plant growth/food production”(Zhao et al. (2021)). Studies have also supported that areas with less anthropogenic activity show low bulk density within the soil. (Cardoso et al. (2013)). Additionally, areas with less human activity tend to have a better accumulation of soil particles, therefore a better structure within the soil. This improved structure is vital for all ecosystems as it allows the permeability for water, air, and roots to be eased, helping them thrive.

There needs to be a shift to sustainable agriculture methods to restore and protect soil health. Sustaining soil health while simultaneously sustaining crop yields is the goal that the agriculture industry should be working towards. In addition, (Katherasala et al. (2024)) emphasize the importance of educating farmers on the harmful effect of current practices and alternative methods. Despite the need for sustainable practices, there is cause for concern and uncertainty towards transitioning to more sustainable agriculture, with most concerns being centered around changes in crop yields and production efficiency. However, (Miner et al. (2020)) provides research showing that under no-tilling conditions, yield declines can be partially mitigated when no-till conditions are combined with residue retention and crop rotation practices. These concerns can also be minimized with the use of soil quality indices (SQIs). SQIs are utilized to measure and monitor soil function and health through chemical and biological properties (Chaudhry et al. (2024)). With practices like crop rotation, residue retention, cover cropping, and livestock integration the SQIs can be determined to monitor soil health with these treatments and ensure it is beneficial for agricultural development. This study looked at different models to determine if crop rotation systems that integrate livestock and cover crops in the Northern Great Plains will result in significantly higher soil quality index scores in relation to microbial activity.

Methods

The data used for our research was collected in the USDA-ARS Northern Great Plains Research Lab from the years 2014-2016. We accessed this data through data.gov from the metadata set called “Conservation Practices Induce Tradeoffs in Soil Function: Observations from the Northern Great Plains” (Liebig et al. (2022)). We downloaded and utilized two of its dataset. The first, SES_On_Station_Intergrative Measures_SMAF, contained different soil scores across different land use types and treatment systems. While the second one, SES_On_Station_Soil_Properties_Near Surface, included measurements for various soil properties for each land use and treatment type. Both of which reflected near-surface (0-5 cm) depth data. After downloading the datasets, we created a new project in RStudio. Once created, we read in each csv file using the read_csv() function, saving each raw dataset to an object. The datasets were then inspected for any shared columns using the glimps() function. We were then able to join the two datasets by their overlapping columns “STUDY”, “LTEXP”, “TRT”, “LAND_USE”, “YEAR”, “REP”, “DEPTH”.

Through Exploratory Data Analysis, we were able to understand the merged data structure. This included identifying outliers, creating visualizations, and getting summary statistics. First we used glimps() to check data types, dimensions, and missing values. To further understand what values were missing. There were two columns that showed a high percentage of missing values, TOC and VESS, to address this we filtered the columns out and created a new dataset with the cleaned data. Then, to view some statistical summaries on our data we used the skim() function, allowing us to view the general trend in mean, standard deviation, and skewness for each variable. As a preprocessor step for correlation analysis, we created a new dataset with only numeric variables. This separation of numeric data allowed us to narrow down the variables that had a strong correlation to the soil quality index score (SMAFQI) variable and to each other. Both vis_cor() and cor() were used for this process. Due to their high correlation with each other, WSA, UREASE, SMAFQI, BGLU, TFAME, PHOS, and BACSUM were all selected to do a separate correlation analysis on. We then visualized each variable’s distribution using the gghistogram() function. Both Beta-Glucosidase Activity (BGLU) and Bacteria Maker FAME (BACSUM) were then selected to complete our remaining research on.

For categorical analysis, we used a box plot, density plot, and violin plots to visualize the relationship between each SMAFQI, BGLU, and BACSUM and treatment (LAND_USE) or management type (TRT). Next, to test for normality in SMAFQI, BGLU, and BACSUM, we used the Shapiro-Wilk normality test. The results informed us that our data did not follow a normal distribution. To help normalize the data, we log-transformed each variable. To verify that the log transformation worked, we used the Shapiro-Wilk normality test again, using the transformed data. Additionally, we viewed the transformation visually using gghistogram(). Both BGLU and BACSUM appeared normally distributed, while SMAFQI still remained skewed. Thus, Pearson’s correlation test was used to test BGLU and BACSUM individually, whereas the Spearman test was used to assess the correlation of the variables without log transforming them. Furthermore, to ensure that SMAFQI, BGLU, and BACSUM were correlated to the different management practices, the Kruskal-Wallis rank sum test was used. To visualize the log-log relationship between BGLU and BACSUM, and view how well both variables predict SMAFQI change, an XY plot of BGLU and BACSUM was made. The scale_color_viridis_c() function was used to color the points by SMAFQI.

We then began a machine learning workflow. First, we set a seed for reproducibility and split the data into an 80% training and 20% testing set. Additionally, a 10-fold cross-validation dataset was created. To preprocess the data recipe() was used to define a series of data preprocessing steps. Next, five different models were built, a workflow_set() was created using the 5 models, and each one was fit to the training data. Their performance was compared using autoplot() and rank_results(). The highest performing model, which was the nearest neighbor model, was then used to build a final workflow on the full training data. The augment() function was then used to make predictions on the test data. To evaluate the final model’s performance metrics() was used and an observed vs predicted scatter plot was made.

Results

# Load and inspect data
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tidymodels)

── Attaching packages ────────────────────────────────────── tidymodels 1.3.0 ──
✔ broom        1.0.8     ✔ rsample      1.3.0
✔ dials        1.4.0     ✔ tune         1.3.0
✔ infer        1.0.8     ✔ workflows    1.2.0
✔ modeldata    1.4.0     ✔ workflowsets 1.1.0
✔ parsnip      1.3.1     ✔ yardstick    1.3.2
✔ recipes      1.3.0     
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ scales::discard() masks purrr::discard()
✖ dplyr::filter()   masks stats::filter()
✖ recipes::fixed()  masks stringr::fixed()
✖ dplyr::lag()      masks stats::lag()
✖ yardstick::spec() masks readr::spec()
✖ recipes::step()   masks stats::step()

integrative_measures_raw <- read_csv("/Users/Avery/github/Final-Project/data/SES_On-Station_Integrative_Measures_SMAF.csv")

Rows: 72 Columns: 21
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): STUDY, LTEXP, TRT, LAND_USE, DEPTH
dbl (16): YEAR, REP, TOCs, AGGs, MBCs, PMNs, PHs, Ps, QCO2s, SBDs, Ecs, AWCs...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

soil_properties_raw <- read_csv("/Users/Avery/github/Final-Project/data/SES_On-Station_Soil_Properties_Near_Surface.csv")

Rows: 72 Columns: 25
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): STUDY, LTEXP, TRT, LAND_USE, DEPTH
dbl (20): YEAR, REP, SBD, POR, SORP, AWHC, WSA, VESS, PH, SLAN, CEC, UREASE,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(integrative_measures_raw)

Rows: 72
Columns: 21
$ STUDY    <chr> "Soil Ecosystem Services", "Soil Ecosystem Services", "Soil E…
$ LTEXP    <chr> "Crop Rotation Study", "Crop Rotation Study", "Crop Rotation …
$ TRT      <chr> "CONSERV", "CONSERV", "CONSERV", "CONSERV", "BAU", "BAU", "BA…
$ LAND_USE <chr> "Continuous cropping", "Continuous cropping", "Continuous cro…
$ YEAR     <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2…
$ REP      <dbl> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1…
$ DEPTH    <chr> "0to5", "0to5", "0to5", "0to5", "0to5", "0to5", "0to5", "0to5…
$ TOCs     <dbl> 0.638, 0.612, 0.718, 0.594, 0.319, 0.398, 0.347, 0.347, 0.634…
$ AGGs     <dbl> 0.857, 0.987, 1.003, 0.977, 0.820, 0.896, 0.918, 0.868, 1.006…
$ MBCs     <dbl> 0.903, 0.924, 0.962, 0.938, 0.986, 0.983, 0.967, 0.976, 1.000…
$ PMNs     <dbl> 0.301, 0.223, 0.478, 0.431, 0.207, 0.217, 0.502, 0.494, 0.763…
$ PHs      <dbl> 0.723, 0.723, 0.723, 0.830, 0.935, 0.918, 0.951, 0.951, 0.976…
$ Ps       <dbl> 1.000, 1.000, 1.000, 0.949, 0.947, 0.962, 0.976, 1.000, 1.000…
$ QCO2s    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ SBDs     <dbl> 0.994, 0.994, 0.994, 0.992, 0.930, 0.988, 0.977, 0.947, 0.919…
$ Ecs      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ AWCs     <dbl> 0.812, 0.830, 0.826, 0.779, 0.731, 0.754, 0.738, 0.708, 0.707…
$ BGs      <dbl> 0.117, 0.102, 0.159, 0.321, 0.326, 0.201, 0.337, 0.340, 1.000…
$ Ks       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ WFPSs    <dbl> 0.978, 0.978, 0.978, 0.974, 0.973, 0.977, 0.976, 0.975, 0.611…
$ SMAFSQI  <dbl> 79.40, 79.79, 83.40, 82.97, 78.26, 79.19, 82.22, 81.58, 89.36…

glimpse(soil_properties_raw)

Rows: 72
Columns: 25
$ STUDY    <chr> "Soil Ecosystem Services", "Soil Ecosystem Services", "Soil E…
$ LTEXP    <chr> "Crop Rotation Study", "Crop Rotation Study", "Crop Rotation …
$ TRT      <chr> "CONSERV", "CONSERV", "CONSERV", "CONSERV", "BAU", "BAU", "BA…
$ LAND_USE <chr> "Continuous cropping", "Continuous cropping", "Continuous cro…
$ YEAR     <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2…
$ REP      <dbl> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1…
$ DEPTH    <chr> "0to5", "0to5", "0to5", "0to5", "0to5", "0to5", "0to5", "0to5…
$ SBD      <dbl> 1.06, 1.02, 1.01, 1.16, 1.25, 1.18, 1.21, 1.24, 1.26, 1.12, 0…
$ POR      <dbl> 0.60, 0.62, 0.62, 0.56, 0.53, 0.55, 0.54, 0.53, 0.52, 0.58, 0…
$ SORP     <dbl> 0.110, 0.168, 0.144, 0.221, 0.074, 0.075, 0.047, 0.049, 0.116…
$ AWHC     <dbl> 0.189, 0.197, 0.196, 0.176, 0.160, 0.168, 0.162, 0.154, 0.153…
$ WSA      <dbl> 33.750, 47.250, 57.750, 45.500, 31.250, 36.750, 38.630, 34.50…
$ VESS     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ PH       <dbl> 4.90, 4.90, 4.90, 5.30, 5.80, 5.70, 5.90, 5.90, 6.10, 6.20, 5…
$ SLAN     <dbl> 157.90, 170.00, 179.60, 189.50, 110.00, 115.30, 107.50, 131.2…
$ CEC      <dbl> 19.0, 20.4, 20.4, 20.2, 17.3, 17.9, 18.4, 18.6, 17.7, 18.4, 2…
$ UREASE   <dbl> 38.19, 44.78, 44.75, 49.47, 68.44, 55.03, 101.47, 75.53, 183.…
$ PHOS     <dbl> 7.15, 11.48, 22.30, 13.24, 33.28, 19.82, 32.96, 27.57, 140.97…
$ BGLU     <dbl> 156.68, 144.09, 186.87, 264.69, 266.40, 211.08, 270.64, 271.8…
$ TFAME    <dbl> 320.14, 333.38, 408.64, 303.57, 343.98, 461.63, 390.40, 345.0…
$ BACSUM   <dbl> 82.32, 84.07, 106.72, 89.28, 91.41, 101.78, 96.70, 96.72, 146…
$ FUNSUM   <dbl> 59.95, 59.63, 69.67, 39.18, 68.59, 103.98, 81.85, 64.00, 123.…
$ NO3N     <dbl> 1.75, 1.95, 3.32, 1.95, 1.95, 1.17, 0.88, 0.88, 31.20, 13.02,…
$ WFPS     <dbl> 0.41, 0.42, 0.41, 0.47, 0.48, 0.44, 0.45, 0.47, 0.73, 0.68, 0…
$ TOC      <dbl> 15.45, 14.49, 15.98, 16.30, 12.89, 13.49, 12.91, 13.26, 18.33…

Source: Article Notebook

# Merge dataset
library(dplyr)

merged_data <- left_join(integrative_measures_raw, soil_properties_raw, by = c("STUDY", "LTEXP", "TRT", "LAND_USE", "YEAR", "REP", "DEPTH"))

# Understand Data Structure – Check types, dimensions, missing values
glimpse(merged_data)

Rows: 72
Columns: 39
$ STUDY    <chr> "Soil Ecosystem Services", "Soil Ecosystem Services", "Soil E…
$ LTEXP    <chr> "Crop Rotation Study", "Crop Rotation Study", "Crop Rotation …
$ TRT      <chr> "CONSERV", "CONSERV", "CONSERV", "CONSERV", "BAU", "BAU", "BA…
$ LAND_USE <chr> "Continuous cropping", "Continuous cropping", "Continuous cro…
$ YEAR     <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2…
$ REP      <dbl> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1…
$ DEPTH    <chr> "0to5", "0to5", "0to5", "0to5", "0to5", "0to5", "0to5", "0to5…
$ TOCs     <dbl> 0.638, 0.612, 0.718, 0.594, 0.319, 0.398, 0.347, 0.347, 0.634…
$ AGGs     <dbl> 0.857, 0.987, 1.003, 0.977, 0.820, 0.896, 0.918, 0.868, 1.006…
$ MBCs     <dbl> 0.903, 0.924, 0.962, 0.938, 0.986, 0.983, 0.967, 0.976, 1.000…
$ PMNs     <dbl> 0.301, 0.223, 0.478, 0.431, 0.207, 0.217, 0.502, 0.494, 0.763…
$ PHs      <dbl> 0.723, 0.723, 0.723, 0.830, 0.935, 0.918, 0.951, 0.951, 0.976…
$ Ps       <dbl> 1.000, 1.000, 1.000, 0.949, 0.947, 0.962, 0.976, 1.000, 1.000…
$ QCO2s    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ SBDs     <dbl> 0.994, 0.994, 0.994, 0.992, 0.930, 0.988, 0.977, 0.947, 0.919…
$ Ecs      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ AWCs     <dbl> 0.812, 0.830, 0.826, 0.779, 0.731, 0.754, 0.738, 0.708, 0.707…
$ BGs      <dbl> 0.117, 0.102, 0.159, 0.321, 0.326, 0.201, 0.337, 0.340, 1.000…
$ Ks       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ WFPSs    <dbl> 0.978, 0.978, 0.978, 0.974, 0.973, 0.977, 0.976, 0.975, 0.611…
$ SMAFSQI  <dbl> 79.40, 79.79, 83.40, 82.97, 78.26, 79.19, 82.22, 81.58, 89.36…
$ SBD      <dbl> 1.06, 1.02, 1.01, 1.16, 1.25, 1.18, 1.21, 1.24, 1.26, 1.12, 0…
$ POR      <dbl> 0.60, 0.62, 0.62, 0.56, 0.53, 0.55, 0.54, 0.53, 0.52, 0.58, 0…
$ SORP     <dbl> 0.110, 0.168, 0.144, 0.221, 0.074, 0.075, 0.047, 0.049, 0.116…
$ AWHC     <dbl> 0.189, 0.197, 0.196, 0.176, 0.160, 0.168, 0.162, 0.154, 0.153…
$ WSA      <dbl> 33.750, 47.250, 57.750, 45.500, 31.250, 36.750, 38.630, 34.50…
$ VESS     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ PH       <dbl> 4.90, 4.90, 4.90, 5.30, 5.80, 5.70, 5.90, 5.90, 6.10, 6.20, 5…
$ SLAN     <dbl> 157.90, 170.00, 179.60, 189.50, 110.00, 115.30, 107.50, 131.2…
$ CEC      <dbl> 19.0, 20.4, 20.4, 20.2, 17.3, 17.9, 18.4, 18.6, 17.7, 18.4, 2…
$ UREASE   <dbl> 38.19, 44.78, 44.75, 49.47, 68.44, 55.03, 101.47, 75.53, 183.…
$ PHOS     <dbl> 7.15, 11.48, 22.30, 13.24, 33.28, 19.82, 32.96, 27.57, 140.97…
$ BGLU     <dbl> 156.68, 144.09, 186.87, 264.69, 266.40, 211.08, 270.64, 271.8…
$ TFAME    <dbl> 320.14, 333.38, 408.64, 303.57, 343.98, 461.63, 390.40, 345.0…
$ BACSUM   <dbl> 82.32, 84.07, 106.72, 89.28, 91.41, 101.78, 96.70, 96.72, 146…
$ FUNSUM   <dbl> 59.95, 59.63, 69.67, 39.18, 68.59, 103.98, 81.85, 64.00, 123.…
$ NO3N     <dbl> 1.75, 1.95, 3.32, 1.95, 1.95, 1.17, 0.88, 0.88, 31.20, 13.02,…
$ WFPS     <dbl> 0.41, 0.42, 0.41, 0.47, 0.48, 0.44, 0.45, 0.47, 0.73, 0.68, 0…
$ TOC      <dbl> 15.45, 14.49, 15.98, 16.30, 12.89, 13.49, 12.91, 13.26, 18.33…

# Identify Data Issues – Outliers, missing data, inconsistencies
library(visdat)
vis_dat(merged_data)

vis_miss(merged_data)

# Clean data
data_clean <- merged_data |>
  select(-VESS, -TOC)

# Summary Statistics – Mean, median, variance, skewness, kurtosis
summary(data_clean)

    STUDY              LTEXP               TRT              LAND_USE        
 Length:72          Length:72          Length:72          Length:72         
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
      YEAR           REP          DEPTH                TOCs       
 Min.   :2014   Min.   :1.00   Length:72          Min.   :0.3190  
 1st Qu.:2014   1st Qu.:1.75   Class :character   1st Qu.:0.5455  
 Median :2015   Median :2.50   Mode  :character   Median :0.6360  
 Mean   :2015   Mean   :2.50                      Mean   :0.6229  
 3rd Qu.:2016   3rd Qu.:3.25                      3rd Qu.:0.7392  
 Max.   :2016   Max.   :4.00                      Max.   :0.8610  
      AGGs             MBCs             PMNs             PHs        
 Min.   :0.4160   Min.   :0.8780   Min.   :0.2070   Min.   :0.5800  
 1st Qu.:0.8905   1st Qu.:0.9928   1st Qu.:0.6813   1st Qu.:0.8400  
 Median :1.0000   Median :0.9990   Median :0.9375   Median :0.9180  
 Mean   :0.9092   Mean   :0.9892   Mean   :0.8122   Mean   :0.8924  
 3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.9925   3rd Qu.:0.9623  
 Max.   :1.0060   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
       Ps             QCO2s        SBDs             Ecs        
 Min.   :0.9020   Min.   :1   Min.   :0.7250   Min.   :0.8480  
 1st Qu.:0.9700   1st Qu.:1   1st Qu.:0.9750   1st Qu.:1.0000  
 Median :1.0000   Median :1   Median :0.9925   Median :1.0000  
 Mean   :0.9848   Mean   :1   Mean   :0.9723   Mean   :0.9971  
 3rd Qu.:1.0000   3rd Qu.:1   3rd Qu.:0.9940   3rd Qu.:1.0000  
 Max.   :1.0000   Max.   :1   Max.   :0.9940   Max.   :1.0000  
      AWCs             BGs               Ks        WFPSs           SMAFSQI     
 Min.   :0.6490   Min.   :0.0630   Min.   :1   Min.   :0.6110   Min.   :75.35  
 1st Qu.:0.7488   1st Qu.:0.2122   1st Qu.:1   1st Qu.:0.9315   1st Qu.:84.09  
 Median :0.7755   Median :0.4705   Median :1   Median :0.9620   Median :88.86  
 Mean   :0.7758   Mean   :0.5118   Mean   :1   Mean   :0.9399   Mean   :87.75  
 3rd Qu.:0.8110   3rd Qu.:0.8648   3rd Qu.:1   3rd Qu.:0.9740   3rd Qu.:93.15  
 Max.   :0.8560   Max.   :1.0080   Max.   :1   Max.   :0.9790   Max.   :95.82  
      SBD             POR              SORP             AWHC       
 Min.   :0.980   Min.   :0.4900   Min.   :0.0440   Min.   :0.1380  
 1st Qu.:1.090   1st Qu.:0.5400   1st Qu.:0.1087   1st Qu.:0.1660  
 Median :1.155   Median :0.5650   Median :0.1380   Median :0.1755  
 Mean   :1.150   Mean   :0.5658   Mean   :0.1649   Mean   :0.1765  
 3rd Qu.:1.210   3rd Qu.:0.5900   3rd Qu.:0.2135   3rd Qu.:0.1890  
 Max.   :1.350   Max.   :0.6300   Max.   :0.5000   Max.   :0.2100  
      WSA              PH             SLAN            CEC       
 Min.   :12.13   Min.   :4.410   Min.   : 65.5   Min.   :13.20  
 1st Qu.:36.31   1st Qu.:5.340   1st Qu.:132.8   1st Qu.:18.48  
 Median :50.76   Median :5.700   Median :177.8   Median :19.90  
 Mean   :48.28   Mean   :5.668   Mean   :184.0   Mean   :19.77  
 3rd Qu.:61.44   3rd Qu.:5.982   3rd Qu.:224.4   3rd Qu.:20.93  
 Max.   :79.00   Max.   :6.700   Max.   :394.5   Max.   :26.10  
     UREASE            PHOS             BGLU            TFAME      
 Min.   : 11.78   Min.   :  7.15   Min.   : 98.62   Min.   :192.6  
 1st Qu.: 44.77   1st Qu.: 59.17   1st Qu.:216.91   1st Qu.:314.9  
 Median : 85.78   Median :125.06   Median :317.92   Median :390.6  
 Mean   : 90.27   Mean   :149.94   Mean   :353.35   Mean   :421.9  
 3rd Qu.:115.15   3rd Qu.:201.10   3rd Qu.:481.12   3rd Qu.:530.8  
 Max.   :216.25   Max.   :585.56   Max.   :878.97   Max.   :959.5  
     BACSUM           FUNSUM            NO3N             WFPS       
 Min.   : 49.18   Min.   : 29.36   Min.   : 0.880   Min.   :0.2800  
 1st Qu.: 83.86   1st Qu.: 57.05   1st Qu.: 6.513   1st Qu.:0.4700  
 Median :103.28   Median : 69.86   Median :11.660   Median :0.5300  
 Mean   :115.69   Mean   : 75.98   Mean   :15.834   Mean   :0.5253  
 3rd Qu.:145.84   3rd Qu.: 95.68   3rd Qu.:19.545   3rd Qu.:0.5825  
 Max.   :281.81   Max.   :150.45   Max.   :80.540   Max.   :0.7300  

library(skimr)
skim(data_clean)

Data summary

Name	data_clean
Number of rows	72
Number of columns	37
_______________________
Column type frequency:
character	5
numeric	32
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
STUDY	0	1	23	23	0	1	0
LTEXP	0	1	19	31	0	3	0
TRT	0	1	3	7	0	2	0
LAND_USE	0	1	19	33	0	6	0
DEPTH	0	1	4	4	0	1	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
YEAR	0	1	2015.00	0.82	2014.00	2014.00	2015.00	2016.00	2016.00	▇▁▇▁▇
REP	0	1	2.50	1.13	1.00	1.75	2.50	3.25	4.00	▇▇▁▇▇
TOCs	0	1	0.62	0.16	0.32	0.55	0.64	0.74	0.86	▅▂▇▆▇
AGGs	0	1	0.91	0.16	0.42	0.89	1.00	1.00	1.01	▁▁▁▁▇
MBCs	0	1	0.99	0.02	0.88	0.99	1.00	1.00	1.00	▁▁▁▁▇
PMNs	0	1	0.81	0.24	0.21	0.68	0.94	0.99	1.00	▁▂▁▂▇
PHs	0	1	0.89	0.10	0.58	0.84	0.92	0.96	1.00	▁▂▁▃▇
Ps	0	1	0.98	0.02	0.90	0.97	1.00	1.00	1.00	▁▁▁▃▇
QCO2s	0	1	1.00	0.00	1.00	1.00	1.00	1.00	1.00	▁▁▇▁▁
SBDs	0	1	0.97	0.05	0.72	0.98	0.99	0.99	0.99	▁▁▁▁▇
Ecs	0	1	1.00	0.02	0.85	1.00	1.00	1.00	1.00	▁▁▁▁▇
AWCs	0	1	0.78	0.04	0.65	0.75	0.78	0.81	0.86	▁▂▇▇▅
BGs	0	1	0.51	0.34	0.06	0.21	0.47	0.86	1.01	▇▃▂▂▇
Ks	0	1	1.00	0.00	1.00	1.00	1.00	1.00	1.00	▁▁▇▁▁
WFPSs	0	1	0.94	0.06	0.61	0.93	0.96	0.97	0.98	▁▁▁▂▇
SMAFSQI	0	1	87.75	5.91	75.35	84.09	88.86	93.15	95.82	▃▂▅▅▇
SBD	0	1	1.15	0.09	0.98	1.09	1.15	1.21	1.35	▃▆▇▆▁
POR	0	1	0.57	0.03	0.49	0.54	0.56	0.59	0.63	▁▆▇▆▃
SORP	0	1	0.16	0.09	0.04	0.11	0.14	0.21	0.50	▇▆▂▁▁
AWHC	0	1	0.18	0.02	0.14	0.17	0.18	0.19	0.21	▁▅▇▅▃
WSA	0	1	48.28	17.69	12.13	36.31	50.76	61.44	79.00	▃▅▇▇▅
PH	0	1	5.67	0.48	4.41	5.34	5.70	5.98	6.70	▁▃▇▇▂
SLAN	0	1	183.99	62.63	65.50	132.78	177.80	224.44	394.50	▆▇▆▂▁
CEC	0	1	19.77	1.90	13.20	18.48	19.90	20.92	26.10	▁▂▇▃▁
UREASE	0	1	90.27	51.72	11.78	44.77	85.78	115.15	216.25	▇▇▆▃▂
PHOS	0	1	149.94	121.37	7.15	59.17	125.06	201.10	585.56	▇▆▂▁▁
BGLU	0	1	353.35	189.55	98.62	216.92	317.92	481.12	878.97	▇▅▆▁▁
TFAME	0	1	421.90	158.00	192.62	314.88	390.56	530.79	959.53	▇▆▃▁▁
BACSUM	0	1	115.69	43.41	49.18	83.86	103.28	145.84	281.81	▇▇▅▁▁
FUNSUM	0	1	75.98	27.51	29.36	57.05	69.85	95.68	150.45	▃▇▃▃▁
NO3N	0	1	15.83	15.89	0.88	6.51	11.66	19.55	80.54	▇▂▁▁▁
WFPS	0	1	0.53	0.09	0.28	0.47	0.53	0.58	0.73	▁▂▇▇▁

# Separating numeric data into a new dataset
numeric_data <- data_clean |>
  select(where(is.numeric))

vis_cor(numeric_data)

Warning in stats::cor(x, method = cor_method, use = na_action): the standard
deviation is zero

data_clean |>
  select(WSA, UREASE, SMAFSQI, BGLU, TFAME, PHOS, BACSUM) |>
  cor()

              WSA    UREASE   SMAFSQI      BGLU     TFAME      PHOS    BACSUM
WSA     1.0000000 0.6957755 0.7090762 0.4889804 0.4685878 0.2274733 0.4523585
UREASE  0.6957755 1.0000000 0.7301251 0.7460760 0.6538916 0.5404495 0.6452244
SMAFSQI 0.7090762 0.7301251 1.0000000 0.7956104 0.5777223 0.6435123 0.5938352
BGLU    0.4889804 0.7460760 0.7956104 1.0000000 0.5505763 0.6265654 0.5533593
TFAME   0.4685878 0.6538916 0.5777223 0.5505763 1.0000000 0.6101661 0.9877583
PHOS    0.2274733 0.5404495 0.6435123 0.6265654 0.6101661 1.0000000 0.6692343
BACSUM  0.4523585 0.6452244 0.5938352 0.5533593 0.9877583 0.6692343 1.0000000

data_clean |>
  select(SMAFSQI, BGLU, BACSUM) |>
  cor()

          SMAFSQI      BGLU    BACSUM
SMAFSQI 1.0000000 0.7956104 0.5938352
BGLU    0.7956104 1.0000000 0.5533593
BACSUM  0.5938352 0.5533593 1.0000000

library(ggpubr)

ggarrange(gghistogram(data_clean$WSA, title = "WSA", bins = 20), gghistogram(data_clean$UREASE, title = "UREASE", bins = 20), gghistogram(data_clean$SMAFSQI, title = "SMAFSQI", bins = 20), gghistogram(data_clean$BGLU, title = "BGLU", bins = 20), gghistogram(data_clean$TFAME, title = "TFAME", bins = 20), gghistogram(data_clean$PHOS, title = "PHOS", bins = 20), gghistogram(data_clean$BACSUM, title = "BACSUM", bins = 20))

# boxplot plot of soil quality index by long-term experiment 
ggplot(data_clean, aes(x = LTEXP, y = SMAFSQI)) +
  geom_boxplot(fill = "skyblue", color = "black") +
  theme_minimal() +
  labs(
    title = "Soil Quality Index by Management Practice",
    x = "Management Practice (LTEXP)",
    y = "SMAFSQI"
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# density plot of variables by management practice/treatment type
ggplot(data_clean, aes(x = SMAFSQI, fill = LTEXP)) +
  geom_density(alpha = 0.5) +
  theme_minimal() +
  labs(
    title = "Density of SMAFSQI by Management Practice",
    x = "Soil Quality Index (SMAFSQI)",
    fill = "Management Practice"
  )

ggplot(data_clean, aes(x = BGLU, fill = LTEXP)) +
  geom_density(alpha = 0.5) +
  theme_minimal() +
  labs(
    title = "Density of SMAFSQI by Management Practice",
    x = "BGLU",
    fill = "Management Practice"
  )

ggplot(data_clean, aes(x = BACSUM, fill = LTEXP)) +
  geom_density(alpha = 0.5) +
  theme_minimal() +
  labs(
    title = "Density of SMAFSQI by Management Practice",
    x = "BACSUM",
    fill = "Management Practice"
  )

# violin plots
SMAFSQI_type_violin <- ggviolin(data_clean, x = "TRT", y = "SMAFSQI", add = "boxplot", color = "TRT")
BGLU_type_violin <- ggviolin(data_clean, x = "TRT", y = "BGLU", add = "boxplot", color = "TRT")
BACSUM_type_violin <- ggviolin(data_clean, x = "TRT", y = "BACSUM", add = "boxplot", color = "TRT")

SMAFSQI_landuse_violin <- ggviolin(data_clean, x = "LAND_USE", y = "SMAFSQI", add = "boxplot", color = "LAND_USE")
BGLU_landuse_violin <- ggviolin(data_clean, x = "LAND_USE", y = "BGLU", add = "boxplot", color = "LAND_USE")
BACSUM_landuse_violin <- ggviolin(data_clean, x = "LAND_USE", y = "BACSUM", add = "boxplot", color = "LAND_USE")

library(patchwork)
combined_violins <- SMAFSQI_landuse_violin / BGLU_landuse_violin / BACSUM_landuse_violin
ggsave("images/combined_violins.png", plot = combined_violins, width = 18, height = 23)

Violin plots (Figure 1) further showed how BGLU, BACSUM, and SMAFSQI were spread out across LAND_USE. In near-surface soil, grazed and ungrazed continuous cropping showed high median values and narrower distributions for SMAFSQI. Likewise, for BGLU and BACSUM, grazed continuous cropping promoted higher β-glucosidase activity and bacterial activity. Hence, Figure 5 illustrates that continuous cropping, especially when grazed, consistently supports higher soil quality. In contrast, across all three indicators (BGLU, BACSUM, and SMAFSQI), spring wheat-fallow systems displayed the lowest medians with a wider spread. Signifying that fallow-based systems may reduce microbial activity and soil quality.

Figure 1

# Assessing Normality
shapiro.test(data_clean$SMAFSQI)

    Shapiro-Wilk normality test

data:  data_clean$SMAFSQI
W = 0.933, p-value = 0.0008572

shapiro.test(data_clean$BGLU)

    Shapiro-Wilk normality test

data:  data_clean$BGLU
W = 0.9281, p-value = 0.0005017

shapiro.test(data_clean$BACSUM)

    Shapiro-Wilk normality test

data:  data_clean$BACSUM
W = 0.9307, p-value = 0.0006649

# normality test on log transformed variables 
shapiro.test(log(data_clean$SMAFSQI))

    Shapiro-Wilk normality test

data:  log(data_clean$SMAFSQI)
W = 0.92437, p-value = 0.0003373

shapiro.test(log(data_clean$BGLU))

    Shapiro-Wilk normality test

data:  log(data_clean$BGLU)
W = 0.97296, p-value = 0.1231

shapiro.test(log(data_clean$BACSUM))

    Shapiro-Wilk normality test

data:  log(data_clean$BACSUM)
W = 0.98818, p-value = 0.7397

ggarrange(gghistogram(log(data_clean$SMAFSQI), title = "SMAFSQI", bins = 20), gghistogram(log(data_clean$BGLU), title = "BGLU", bins = 20), gghistogram(log(data_clean$BACSUM), title = "BACSUM", bins = 20))

# Correlation

# Pearson’s correlation test, works on normally distributed data
cor.test(log(data_clean$BGLU), log(data_clean$BACSUM))

    Pearson's product-moment correlation

data:  log(data_clean$BGLU) and log(data_clean$BACSUM)
t = 6.2409, df = 70, p-value = 2.928e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4251322 0.7286502
sample estimates:
      cor 
0.5979125 

# Spearman correlation test, works on non-parametric and non-linear data sets
cor.test(data_clean$BGLU, data_clean$BACSUM, method = "spearman")

    Spearman's rank correlation rho

data:  data_clean$BGLU and data_clean$BACSUM
S = 28654, p-value = 1.544e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.5392951 

cor.test(data_clean$BGLU, data_clean$SMAFSQI, method = "spearman")

Warning in cor.test.default(data_clean$BGLU, data_clean$SMAFSQI, method =
"spearman"): Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  data_clean$BGLU and data_clean$SMAFSQI
S = 7232.6, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.8837135 

cor.test(data_clean$SMAFSQI, data_clean$BACSUM, method = "spearman")

Warning in cor.test.default(data_clean$SMAFSQI, data_clean$BACSUM, method =
"spearman"): Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  data_clean$SMAFSQI and data_clean$BACSUM
S = 25712, p-value = 6.137e-08
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.5866019 

Spearman’s rank correlation revealed moderately strong relationships between the variables BGLU, BACSUM, and SMAFSQI. Beta-Glucosidase Activity (BGLU) had a strong positive correlation with Soil Quality Index Score (SMAFSQI) (ρ = 0.88, p < 2.2e-16), suggesting that greater β-glucosidase activity was associated with improved soil quality. Bacteria Marker FAME (BACSUM) had a moderate positive correlation with SMAFSQI (ρ = 0.59, p < 6.14e-08), proposing that increasing bacterial FAME could lead to an increase in soil quality. Similarly, BGLU and BACSUM were moderately correlated with each other (ρ = 0.54, p < 1.54e-06). Each relation had a p-value of < 0.05, indicating that each correlation was statistically significant (Figure 2).

Figure 2

# Kruskal-Wallis rank sum test
kruskal.test(SMAFSQI ~ LAND_USE, data = data_clean)

    Kruskal-Wallis rank sum test

data:  SMAFSQI by LAND_USE
Kruskal-Wallis chi-squared = 56.642, df = 5, p-value = 5.993e-11

kruskal.test(BGLU ~ LAND_USE, data = data_clean)

    Kruskal-Wallis rank sum test

data:  BGLU by LAND_USE
Kruskal-Wallis chi-squared = 43.918, df = 5, p-value = 2.407e-08

kruskal.test(BACSUM ~ LAND_USE, data = data_clean)

    Kruskal-Wallis rank sum test

data:  BACSUM by LAND_USE
Kruskal-Wallis chi-squared = 35.721, df = 5, p-value = 1.08e-06

# Kruskal-Wallis rank sum test
kruskal.test(SMAFSQI ~ LTEXP, data = data_clean)

    Kruskal-Wallis rank sum test

data:  SMAFSQI by LTEXP
Kruskal-Wallis chi-squared = 52.656, df = 2, p-value = 3.68e-12

kruskal.test(BGLU ~ LTEXP, data = data_clean)

    Kruskal-Wallis rank sum test

data:  BGLU by LTEXP
Kruskal-Wallis chi-squared = 41.85, df = 2, p-value = 8.175e-10

kruskal.test(BACSUM ~ LTEXP, data = data_clean)

    Kruskal-Wallis rank sum test

data:  BACSUM by LTEXP
Kruskal-Wallis chi-squared = 25.8, df = 2, p-value = 2.498e-06

Kruskal-Wallis tests indicated that BGLU, BACSUM, and SMAFSQI significantly varied across different land use practices (LAND_USE) and management treatments (LTEXP). Across land use practices, SMAFSQI differed significantly for each practice (p = 5.99e-11), along with BGLU (p = 2.41e-08) and BACSUM (p = 1.08e-06) (Figure 3). Similarly, significant variability was found across management treatments for SMAFSQI (p = 3.68e-12), BGLU (p = 8.18e-10), and BACSUM (p = 2.50e-06) (Figure 3). The test results support the assumption that management strategies affect both soil properties and soil quality.

Figure 3

log_visualization <- ggplot(data_clean, aes(x = BGLU, y = BACSUM)) +
  geom_point(aes(color = SMAFSQI)) +
  geom_smooth(method = "lm") +
  scale_color_viridis_c() +
  # Apply log transformations to the x and y axes
  scale_x_log10() + 
  scale_y_log10() +
  theme_linedraw() +
  theme(legend.position = "bottom") + 
  labs(title = "Beta-Glucosidase Activity vs Bacteria-marker FAME vs Soil Quality Index Score", 
       x = "Beta-Glucosidase Activity", 
       y = "Bacteria-marker FAME",
       color = "Soil Quality Index Score")

ggsave("images/log_visualization.png", plot = log_visualization, width = 8, height = 7, units = "in")

`geom_smooth()` using formula = 'y ~ x'

The scatterplot of BGLU and BACSUM, whose points were colored by SMAFSQI, visually represented these correlations. The plot showed a gradient of soil quality scores. The higher soil quality scores generally clustered where both BGLU and BACSUM were high. Whereas the lower soil quality scores were gathered in areas of low BGLU and BACSUM (Figure 4). Signifying that BGLU and BACSUM could be used to predict the soil quality score.

Figure 4

set.seed(123)

data_split <- initial_split(data_clean, prop = 0.8)
train_data <- training(data_split)
test_data  <- testing(data_split)

cv <- vfold_cv(train_data, v = 10)

soil_recipe <- recipe(SMAFSQI ~ BGLU + BACSUM, data = train_data) |>
  step_log(all_predictors()) |>
  step_interact(terms = ~ BGLU:BACSUM) |>
  step_naomit(all_predictors(), all_outcomes())

lm_model <- linear_reg() |>
  set_engine('lm') |>
  set_mode("regression")

rf_model <- rand_forest() |>
  set_engine("ranger", importance = "impurity") |>
  set_mode("regression")

bt_model <- boost_tree() |>
  set_engine("xgboost") |>
  set_mode("regression")

bag_model <- bag_mlp() |>
  set_engine("nnet") |>
  set_mode("regression")

nn_model <- nearest_neighbor() |>
  set_engine("kknn") |>
  set_mode("regression")

library(tidyverse)
library(tidymodels)
library(powerjoin)
library(glue)
library(vip)

Attaching package: 'vip'

The following object is masked from 'package:utils':

    vi

library(baguette)

wf_set <- workflow_set(list(soil_recipe), list(lm_model, rf_model, bt_model, bag_model, nn_model)) |>
  workflow_map('fit_resamples', resamples = cv)

autoplot(wf_set)

rank_results(wf_set, rank_metric = "rsq", select_best = TRUE)

# A tibble: 10 × 9
   wflow_id         .config .metric  mean std_err     n preprocessor model  rank
   <chr>            <chr>   <chr>   <dbl>   <dbl> <int> <chr>        <chr> <int>
 1 recipe_nearest_… Prepro… rmse    1.97   0.304     10 recipe       near…     1
 2 recipe_nearest_… Prepro… rsq     0.913  0.0277    10 recipe       near…     1
 3 recipe_bag_mlp   Prepro… rmse    3.17   0.735     10 recipe       bag_…     2
 4 recipe_bag_mlp   Prepro… rsq     0.876  0.0279    10 recipe       bag_…     2
 5 recipe_rand_for… Prepro… rmse    2.43   0.304     10 recipe       rand…     3
 6 recipe_rand_for… Prepro… rsq     0.862  0.0377    10 recipe       rand…     3
 7 recipe_linear_r… Prepro… rmse    2.37   0.381     10 recipe       line…     4
 8 recipe_linear_r… Prepro… rsq     0.859  0.0431    10 recipe       line…     4
 9 recipe_boost_tr… Prepro… rmse    3.03   0.327     10 recipe       boos…     5
10 recipe_boost_tr… Prepro… rsq     0.824  0.0335    10 recipe       boos…     5

final_fit <- workflow() |>
  add_recipe(soil_recipe) |>
  add_model(nn_model) |>
  fit(data = train_data)

final_preds <- augment(final_fit, new_data = test_data)

metrics(final_preds, truth = SMAFSQI, estimate = .pred)

# A tibble: 3 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.98 
2 rsq     standard       0.931
3 mae     standard       1.68 

# Create a plot of the observed vs predicted values with clear title, axis labels, and a compelling color scale
final_model_fit <- ggplot(final_preds, aes(x = SMAFSQI, y = .pred, colour = SMAFSQI)) +
  scale_color_viridis_c() +
  geom_point() +
  geom_abline() +
  theme_linedraw() +
  labs(title = "Observed vs Predicted Soil Quality Index Score",
       x = "Observed",
       y = "Predicted")

ggsave("images/final_model_fit.png", plot = final_model_fit, width = 10, height = 7, units = "in")

Among the five machine learning models evaluated, the nearest neighbor model showed the best performance. The final fitted model had an RMSE of 1.98 and an R² of 0.93. These results suggest a strong relationship between BGLU and BACSUM as predictors for SMAFSQI. The observed vs. predicted scatterplot (Figure 5) exhibited a clustering of data points along the 1:1 reference line, thus indicating a very strong model fit with minimal deviation.

Figure 5

Discussion/Conclusion

Our main findings support the hypothesis that crop-livestock integration and the use of cover crops result in significantly higher soil quality index scores in rainfed Northern Great Plains systems. These conservation practices stimulate microbial activity and enhance soil biological structure, offering a resilient path forward for sustainable dryland agriculture. Diversification and cover cropping support soil function and microbial activity. With these conservation practices, the increase of BGLU and BACSUM is observed, leading to the conclusion that the two variables can be used to predict the soil quality score. Additional findings from our research tell us that fallow-based methods are linked to a decrease in microbial activity. This is a method that is widely utilized in agricultural areas, as it restores organic matter, but this restoration ultimately falls below the performance of continuous cropping and diversification.

Other research from (Miner et al. (2020)) suggests that environmental benefits from crop diversity are regionally specific, with agroecosystems globally being very different. Due to this, it cannot be said for all agricultural systems that diversification and the use of cover crops will increase yield production and SQIs. However, for the rainfed croplands of Bismarck, North Dakota, this is the expected outcome. Our data and study suggest that biological pathways, like microbial stimulation from plant diversity, could support yield improvements. This microbial stimulation from plant diversity is supported by others, as Jalloh states, as it is associated with decomposition, carbon utilization, and nitrogen fixation, all key components of helping the stability of an ecosystem (Jalloh et al. (2024)).

The major limitation of this study was the small sample size of the data we had to work with. There were many difficulties at the beginning of our study that led us to change the scope of the study itself. Due to this small sample size, we were unable to see the desired results through our preliminary testing. Many tests fell short, specifically correlation tests, with there simply not being enough data in the originally used datasets. Because of this issue, we were forced to change our original hypothesis to align with newly found datasets. Another issue we encountered with the datasets was the time scale they offered. While one dataset we used had sufficient data for multiple years, other datasets with very similar data would have inputs from only one year. This lack of consistency proved to be more of a challenge for us than we anticipated.

For future studies, we recommend a larger sample size. As well as expanding the study outside of the small cropland in North Dakota. The smaller spatial scale was beneficial for our study, as we had limited time to collect data and work with it, but to determine which conservation practices would be beneficial for other types of biomes and ecosystems, more studies specific to different locations need to be done, as one method that works for one ecosystem may not work for another.

This study supports shifting from traditional practices to conservation practices like cover cropping and diversified rotations that maintain or improve SQIs for long-term sustainability. Although concerns for crop yield and production are widespread, yield efficiency is likely to increase as the soil quality and health improve. The SQIs can be utilized to determine what will help soil health and function in an agricultural system, which will, in turn, help the production and efficiency of crops.

References

Belete, T., & Yadete, E. (2023). Effect of mono cropping on soil health and fertility management for sustainable agriculture practices: A review. Plant Sci, 11, 192–197.

Cardoso, E. J. B. N., Vasconcellos, R. L. F., Bini, D., Miyauchi, M. Y. H., Santos, C. A. dos, Alves, P. R. L., Paula, A. M. de, Nakatani, A. S., Pereira, J. de M., & Nogueira, M. A. (2013). Soil health: Looking for suitable indicators. What should be considered to assess the effects of use and management on soil health? Scientia Agricola, 70, 274–289.

Chaudhry, H., Vasava, H. B., Chen, S., Saurette, D., Beri, A., Gillespie, A., & Biswas, A. (2024). Evaluating the soil quality index using three methods to assess soil fertility. Sensors, 24(3), 864.

Jalloh, A. A., Khamis, F. M., Yusuf, A. A., Subramanian, S., & Mutyambai, D. M. (2024). Long-term push–pull cropping system shifts soil and maize-root microbiome diversity paving way to resilient farming system. BMC Microbiology, 24(1), 92.

Katherasala, S., Bheenaveni, R. S., Chinthakindi, P., Vadlakonda, D., Kummari, N., & Bandi, S. S. (2024). Unveiling the detrimental impacts of intensive chemical use and monoculture on soil health and sustainable development goal 15: Life on land in telangana. Journal of Lifestyle and SDGs Review, 4(2), 10–47172.

Liebig, M. A., Acosta-Martinez, V., Archer, D. W., Halvorson, J. J., Hendrickson, J. R., Kronberg, S. L., Samson-Liebig, S. E., & Vetter, J. M. (2022). Conservation practices induce tradeoffs in soil function: Observations from the northern great plains. Soil Science Society of America Journal, 86(6), 1413–1430.

Miner, G. L., Delgado, J. A., Ippolito, J. A., & Stewart, C. E. (2020). Soil health management practices and crop productivity. Agricultural & Environmental Letters, 5(1), e20023.

Worley, N. (2024). The truth of industrial agriculture. The Chanterelle, 1(1), 9.

Zhao, F., Wu, Y., Hui, J., Sivakumar, B., Meng, X., & Liu, S. (2021). Projected soil organic carbon loss in response to climate warming and soil water content in a loess watershed. Carbon Balance and Management, 16, 1–14.